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Abstract 

We study the long-range effective drift and diffusivity of a particle in a 
random medium moving subject to a given molecular diffusivity and a local 
drift. The local drift models the effect of a random electrostatic field on a 
neutral but polarizable molecule. Although the electrostatic field is assumed 
to obey Gaussian statistics the induced statistics of the drift velocity field are 
non-Gaussian. 

We show that a four-loop perturbation theory calculation of the effective 
diffusivity is in rather good agreement with the outcome of a numerical sim- 
ulation for a reasonable range of the disorder parameter. We also measure 
the effective drift in our simulation and confirm the validity of the "Einstein 
relation" that expresses the equality of the renormalization factors, induced 
by the random medium, for the effective drift and effective diffusivity, rela- 
tive to their molecular values. The Einstein relation has previously only been 
confirmed for Gaussian random drift fields. The simulation result, for our 
non-Gaussian drift model, is consistent with a previous theoretical analysis 
showing the Einstein relation should remain true, independently of the precise 
character of the statistics of the drift velocity field. 
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1 Introduction 



The advective diffusion of passive scalar fields in random environments has been 
extensively studied by both analytical and numerical techniques, with particular 
emphasis on the computation of effective parameters for a diffusion process that 
combines molecular diffusion with a drift term that depends linearly on the gradient 
of a random scalar field |I|| - pOU - 

The problem is well studied and understood in the case of transport by a gradient 
velocity field which exhibits Gaussian statistics. In isotropic systems a renormaliza- 
tion group approach can be shown to give exact results in one and two dimensions 
[|I], ^ 0. Somewhat surprisingly the same approach works extremely well for the 
isotropic problem in three dimensions ||, 0, though the situation is less clear 
in the absence of isotropy P, |[. 

One ingredient in the success of the renormalization group method is that, at 
each stage of the calculation, it respects the Einstein relation that guarantees the 
equality of the renormalization factors of the effective diffusivity, k^, and drift, Ag, 
parameters relative to their molecular values, kq and Aq repectively. That is 

Ag 

Kq Aq 

This result does not hold in the, somewhat arbitrary, Hartree-Fock resummation 
procedure which turns out to be even less accurate than simple low order perturba- 
tion theory which also breaks down for strong disorder. 

It turns out that the Einstein relation holds in a wide range of circumstances 
including ones where the diffusivity and drift coefficient have non-trivial tensorial 



structure |T^. The only requirement is that the two tensors are linearly depen- 
dent. Thus if at the molecular level we have Aq ij = tkq ij then we will find for the 
macroscopic effective parameters that Ag ij = tk^ ij . The isotropic situation is an 
example of this. The reason that the proportionality factor survives renormalization 
can be traced back to the existence of a finite sample equilibrium distribution with a 
vanishing micro- current. In equilibrium statistical mechanics this vanishing is guar- 
anteed which explains the emergence of the Einstein relation in this context. In 
the quenched models with which we are concerned the Einstein relation is not guar- 
anteed except in the special but important circumstances indicated above. Given 
these conditions of proportionality of the molecular tensors however, the Einstein 
relation can be shown to hold sample by sample. It follows that the relation will hold 
independent of the precise nature of the statistics of the gradient velocity field. This 
theoretical prediction lends considerable interest to the results of the numerical 



simulation of the physically motivated model that we study in this paper, in which 
the gradient velocity field does not have Gaussian statistics. 

A physical realization of the Gaussian case is the diffusion of an ion in a sea of 
fixed, disordered charges that give rise to an electrostatic potential. Such a situation 
can be created by embedding fixed ions in the rather accommodating structure of a 
zeolite matrix. It is not unreasonable to treat the resulting electrostatic field as hav- 
ing Gaussian statistics at length scales somewhat larger than the molecular level but 
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very short compared to macroscopic scales. A related problem, particularly relevant 
to applications in chemical processes, concerns the diffusion of a neutral molecule 
in such a medium. The case of the diffusion of benzene in zeolites with random 
ionic inclusions has been studied experimentally and has important technological 
applications |Tl|. In this situation the force on the molecule arises because of the 



interaction of its induced dipole moment with the electrostatic field produced by the 
included ions. The result is a gradient velocity field with non-Gaussian statistics. 
For simplicity we take it for granted that all tensorial effects are absent. 

We assume that the dipole moment P, of the diffusing particle is proportional to 
the ambient electric field Ei thus 

Pi = fiEi . (2) 

The force on the molecule is 

F, = Pjd,E, , (3) 
and since for an electrostatic field djEi = diEj , it follows that 

F, = fiEjdiEj = ^lidiE] . (4) 

If we introduce the electrostatic potential </'(x) then the force on the molecule be- 
comes 

= ^M(V0(x))2 _ 
If we assume that the drift of the molecule is given by 

Mi(x) = uFi , (6) 
then, absorbing all constants into an overall parameter Aq, we have 

M,(x) = ^Ao9,(V0(x))2 . (7) 

It is still reasonable to model the statistical properties of the electrostatic potential 
0(x) by a Gaussian field. However because of the quadratic relationship between 
Mj(x) and 0(x) the statistics of the drift velocity are not Gaussian in this model. 
The higher cumulants of the velocity field beyond the second do not vanish. Indeed 
it is very clear that the third and other odd order cumulants exist. This means that 
the problem is not symmetrical under the change Xq —>■ —Xq . For this reason we 
explore both positive and negative values of Aq even though the physical derivation 
of the model suggests that Aq > . For the purposes of the simulation we take the 
opportunity to absorb the normalization of 0(x) into the definition of Aq and require 

m^))') = 1 • (8) 

The problem we address then, is the evaluation of the long-range effective pa- 
rameters in terms of the (non-Gaussian) statistical properties of the random flow. 
Such a model presents its own new technical difficulties. Calculational schemes such 
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as Self-Consistent Perturbation Theory or the Renormahzation Group (RG), which 
proved so successful in Gaussian problem, turn out to be hard to apply to this 
particular problem. The barriers to a straightforward application of these more so- 
phisticated perturbation methods arise mainly because of the increased complexity 
of the vertex structure at low wave number that is quickly revealed by perturbation 
theory. In addition the lowest order correction to the propagator begins at two loops 
rather than one loop as in the Gaussian case. In this paper we confine ourselves 
to computing the standard perturbation expansion for the effective diffusivity to 
four-loop order. As will become clear the results are consistent with the outcome of 
our numerical simulation of the model over a significant range of values for the drift 
parameter. 

2 Perturbation Theory 

The diffusion equation in which we are interested has the form 

^^!^ = V(«:oVP(x,t)-u(x)P(x,t)) . (9) 

Here, P(x, t) is the probability density of a particle moving according to the equation 

x = u(x)+w(t) , (10) 

where w(t) is a white noise term that satisfies 

{wi{t)) = and {wi{t)wj{t')) = 2Ko6ij6{t - t') . (11) 

The flow field, u(x), is taken to be time independent and is the gradient of a scalar 
field 

u(x) = AoV^(x) , (12) 

but, because the flow originates in the interaction of the induced dipole moment 
of the diffusing particle with the electrostatic field, ip{x.) does not exhibit Gaussian 
statistics. In fact, as indicated above, we have 

^(x) = i(V</>(x))2 , (13) 

where 0(x) is a homogeneous Gaussian random field characterized by the disorder 
averages 

((/.(x)) = and (0(x)(/.(y))=A(x-y) . (14) 

For simplicity we take the disorder to be isotropic, that is, A = A(|x|). 

The perturbative approach to solving equation (1) is well known ||l], ||, |^, |, 0, 
1^, |15| and we only summarize here the necessary results. Since we are interested 
in the effective parameters governing the evolution of the distribution P(x, t), we 
study the related static Green function, G'(x), which satifies 

KoV'G(x) - V(u(x)G(x)) = -5(x) . (15) 
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A perturbation series in the coupling Aq for G(k) can be generated by iterating the 
formal solution to equation (7) in Fourier space: 

The Green function averaged over the velocity ensemble, (G(k)), can be written as 

where the averaging over the velocity ensemble is done using Wick's theorem to give 
a diagrammatic expansion and S(k) is the summation of one particle- irreducible 
diagrams. The expected asymptotic behaviour of the diffusion process at large 
distances and times implies that the small k behaviour of (G(k)) is given by 

ft^e = fi^o - ^S(k)|k=o , (18) 

where Kg is the effective diffusivity. The Feynman rules for the diagrammatic ex- 
pansion are as follows: 

1. Wave vector is conserved at each vertex; 

2. Each full line carries a factor l/^ok^; 

3. Wave vector is integrated around closed loops with a factor c?^q/(27r)^; 

4. Each vertex, whose diagrammatic representation is shown in figure |l|, carries 
a factor — Aq q ■ p k ■ (q + p) ; 

5. Each internal dashed line carries a factor A(q); 

6. Each diagram must be divided by the usual symmetry factor. 
In what follows, we use the exphcit spectrum 

A(q) = ^^e-'''/''« • (19) 
fco 

The normalization is chosen so that ((0(x))^) = 1. In our numerical calculations, 
we set fco = 1- 

There is no one loop correction to the propagator. As mentioned before, this 
together with the fact that new vertices are generated when correcting the primitive 
one makes it extremely delicate to implement other perturbative schemes such as 
Self-Consistent or RG methods. Therefore, we concentrate on a straightfoward 
perturbation theory calculation which has to be carried out to at least four-loop 
order to get a sensible outcome. The formal manipulations of which we make use 
are basically the same as those utilized in [§, ^ |^ and we only elucidate the more 
elaborate steps and state the main results. 
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The two-loop contribution to S(k) is associated with the diagram in figure 0. 
According to the above Feynman rules it is 

E(^)(k) = -i^ / 7^7^ A(q) A(p) ■ q)' • (p + q) ('^ - V p) • + p) , 

^ ' 2kJ (27r)3(27r)3 (k-q-p)^ 

(20) 

and it can be easily computed to O(k^) with the result 

E(2)(k) = l^k2 . (21) 
2 kq 

In order to encounter deviations from Gaussian behaviour we must calculate 
beyond second order and consider the three-point correlator of the velocity field, 
that is, to include the first non-zero odd power in Aq in the perturbation expansion. 
The three-loop diagram contributing to S(k) is shown in figure |] and yields 

Aq f d^q cPp d^t 
4 J (27r)3 (27r)3 (27r)3 
k-(q + p)(k-q-p)-(t-p)(k-q-t)-(q + t) 



(k-q-p)2(k-q-tf 

The numerator can be re-written using the identity 
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[k-q- t) ■ (-q- t) = (k-q- t)^ -k- (k-q- t) . (23) 



We then have 



k - (q + p) (k-q-p) ■ (p-1 
k-(k-q-t) 1 



(k-q-p)^ 



(k-q-p)' (k-q-t)^ 



(24) 



The first term, which turns out to be the dominant one, can be performed analyti- 
cally whereas the second one, much smaller, has to be done numerically to O(k^). It 
is instructive to explain in some detail how to compute the analytical contribution 
to S*^'^)(k). In the first term of equation (pij), the term odd in t integrates to zero 
and, moreover, it can be symmetrized for p and q to yield the result 

IXl f d^q d^p dh 
24J (27r)3 (27r)2 {2n) 
k ■ (q + p) (k - q - p) ■ (q + p) 



(25) 



(k-q-p 

Using again equation ( ^31) (for vectors p and q) leads to 



4 J (27r)3 (27r)3 (27r)3 

- q- : 
q- p) 
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The first term in equation ( pq) integrates to zero. We need to evaluate the second 
term only to O(k^). Because of the explicit factors of k in the integrand we can set 
k = everywhere else to obtain 



c/^q d^p dH 
(27r)3 (27r)3 (27r)3 
k- (q + p)k- (q + p) 



1A| 
2 Kq 



A(q)A(p)A(t)qpptqt 



[q + p) 



(27) 



This is easily evaluated as 



Siink) = i^k^ 



2 Kq 



d' 



q 



1 3 



(2vr) 



q'A(q) 



2 Kg 



Combining the analytical and numerical pieces yields the following result 

S(3)(k) = ^k^ (1 + 0.030375) . 



(2^ 



(29) 



This last contribution is the first of the odd power terms in the expansion that 
are responsible for the asymmetry of Kg under change of sign of the coupling Aq . Its 
presence is a direct confirmation of the non-Gaussian property of the statistics of the 
velocity field. However, truncating the power series at 0{Xl) causes Kp to increase 
for larger negative values of Aq, as is shown in figure ^. This unphysical feature is 
obviously an artifact of perturbation theory and can be circumvented by including 
the fouth-order term in the perturbation expansion. The diagrams corresponding 
to this order are shown in figure 11 and give the following contributions: 



E(^'^)(k) 



lA^ f d\ d^q' d^p d^p' ~ 



A(q)A(q')A(p)A(p') (q-q')'(p-p' 



4 Kg 7 (27r)3 (27r)3 (27r)3 (27r)3 



k-(p + p') (k-p-p')-(q + q') 

_ (k - p - p' - q - qQ • (p + pQ (k - q - qQ ■ (q + q') 
(k - p - p')^ (k - p - p' - q - q')^ (k - q - q')^ 



(30) 



E(^'')(k) 



lA^ r d\ d^q' d^p d^p' r 



A(q)A(q')A(p)A(p') (q-q')'(p-p' 



44 J (27r)3 (27r)3 (27r)3 (27r)3 
k-(p + p') (k-p-p')-(q + q') 
(k - p - p^ - q - qQ ■ (q + qQ (k - p - p') ■ {p + p' 
(k - p - p')^ (k - p - p' - q - q')^ 



(31) 



S(^")(k) 



A^ 



d^q d^q' d^p d^p' - 



kI J (27r)3 (27r)3 (27r)3 (27r)3 
p ■ qp - q'k- (p + p) (k - p- p') ■ (q- p') 
(k - p - q) ■ (q^ - q) (k - p - qQ ■ (p + q^ 
(k - p - p')^ (k - p - q)^ (k - p - q')^ 



A(q)A(q')A(p)A(pOq-q'p-p' 



(32) 
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Kj] J (27r)3 (27r)3 (27r)3 (27r)3 

p' ■ q p ■ q' k ■ (p + p') (k - p - pO ■ (q - p') 
_ (k - p - q) ■ (q^ - p) (k - q - q') ■ (q + qQ 
(k - p - p')^ (k - p - q)^ (k - q - q')^ 

f d^q d^q' c/3p ^3p/ 



(33) 



€f] y (27r)3 (27r)3 (27r)3 (27r)3 
p' ■ q'p • qk - (p + p) (k - p - p') • (q+ q') 
_ (k - q - q^ - p - pQ ■ (q^ + p') (k - q - p) ■ (q + p) 
(k - p - p')^ (k - q - q' - p - p')^ (k - q - p)^ 



(34) 



The above expressions can be simplified using the same sort of manipulations as 
in the three-loop case. The first two are calculated analytically to O(k^) with the 
result 

'1 f rf3q 

The remaining ones lead, again, to a dominant contribution that can be calculated 
analytically plus smaller pieces which are evaluated numerically to yield 



Kn 



1\' 



^{A,c+d+e) = ^ (- - 0.06) . (36) 

Kq 2 



The outcome for Kg to 0{Xq) is then 



fi:e = «:ojl- ^^- 0.530375^ -0.315^1 . (37) 

I 2 Kg Kg \ 



The results for Kg at two, three and four loops are shown plotted in figure ^ 
for the range —1.5 < Aq < 1.5 and kq = 1. Clearly, the four-loop perturbative 
calculation is successful in surmounting the difficulties of the three-loop case while 
encoding the deviations from Gaussian behaviour. The latter, are translated into a 
fast decay of Kg for positive values of Aq whereas it decreases more slowly for the 
negative values. 



3 Numerical Simulation of Drift and Diffusivity 

To simulate the evolution of the scalar field P(x, t) we integrate numerically the 
stochastic equation for the evolution of a particle with path x(t) given by equation 
([T0[) . The resulting probability distribution for the particle position x(t) is then 
P(x, t) with the initial condition P(x, 0) = 5(x). 

The discrete form of equation ([T0|) suitable for numerical integration is 

x„,+i -x„ = u(x„)At+ (2KoAt)i/2^„ , (38) 
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where ^„ is a Gaussian random three- vector of zero mean and unit variance for each 
component. This equation models equation ([T0|) correctly to 0{At) but the details 
of a third-order Runge-Kutta scheme correct to 0{At^) are given in fl^. We use 
this third-order scheme in our numerical simulation. 

The realizations of the random field 0(x) are constructed in the usual way 
g,[T|. We set 

2x1/2 N 

— j ^ cos(k„ • X + e„) , (39) 

' n=l 

where the vector e„ is distributed uniformly over the unit sphere and the wavevector 
k„ is distributed according to the distribution 

For sufficiently large the central limit theorem guarantees that 0(x) is Gaussian 
up to 0{1/N) corrections. 

The effective diffusivity is computed, for a realization of the velocity field, from 
the ensemble of paths by 

1 M 



a=l 



= 6M + 0(1) as t->cx) . (41) 

To measure the effective drift, Ag, we add a constant drift term to equation (4). 
In appropriate units, it is given by 

u' = Aog , (42) 

where g is a uniform gradient field. Assuming that the latter lies along the x-axis, 
for a realization of the velocity field, Ag is computed according to 

1 M 

(x(t)) = lim — Vx"(t) 

a=l 

= Xegt + 0{1) as t^oo . (43) 

In practise, the number of field realizations and M are finite but large enough 
to give an estimate of Kg and Ag with reasonable error. In addition the simulation 
must be carried to values of t large enough to ensure that measurements are being 
performed in the asymptotic regime controlled by the long-range or "renormalized" 
parameters. This is tested by ensuring that the estimates for Kg and Ae are indepen- 
dent of the range of t used evaluate them, within statistical errors. In our simulation 
we tracked the trajectories of 400 particles in each of 1280 realizations of the velocity 
field, each of them containing 128 modes. The simulation was run for a total of 8000 
(32000) time steps of length 0.025 (0.0125) for the smaller (larger) absolute values 
of the bare coupling Aq. For larger values of Aq and because of time limitations the 
number of paths followed was reduced to a minimum of a hundred. 
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Ao 


g 




(Ae/Ao) 


-1.5 


0.05 


0.69276(60) 


0.6946(9) 


-1.2 


0.05 


0.75797(66) 


0.7583(12) 


-1.0 


0.05 


0.80412(69) 


0.8032(15) 


-0.7 


0.05 


0.87360(75) 


0.8753(31) 


-0.5 


0.05 


0.92560(56) 


0.9213(32) 


-0.2 


0.125 


0.98386(59) 


0.9803(33) 


0.2 


0.125 


0.97632(59) 


0.9790(33) 


0.5 


0.05 


0.79788(49) 


0.7983(30) 


0.7 


0.05 


0.54744(52) 


0.5458(25) 


1.0 


0.05 


0.19537(25) 


0.1939(8) 


1.2 


0.05 


0.08979(15) 


0.0895(5) 


1.5 


0.05 


0.03059(7) 


0.0296(3) 



Table 1: Measurements of Kg and Ag for various values of the disorder parameter Aq. 



The results from the perturbation theory calculation of against results from 
the simulation are shown in figure |^ for the range —1.5 < Aq < 1.5 and assuming 
Kq = 1. The simulation confirms the pronounced asymmetry of Kg as a function of Aq, 
exhibiting a fast (slow) decay for positive (negative) values of Aq. The four-loop order 
results compare well with the simulation outcome in the range —0.7 < Aq < 0.7. 
As expected, the four-loop perturbative calculation encodes the relevant qualitative 
features of the flow and is effective in the region of small disorder as opposed to the 
two-loop and three-loop case. 

In table 1 we show the measurements of both Ke and Ae over a range of values of 
the disorder parameter (again, we take kq = 1). The results clearly show that the 
equality of the two renormalized factors is well maintained throughout with only 
slight discrepancies some cases due to systematic errors. For the higher values of 
the disorder parameter another possible source of error is that the value of the drift 
parameter has become so large that 0{g^) effects are influencing the values of the 
measured quantities. Nevertheless, we are confident that the simulation supports the 
conclusion that in gradient flow, irrespective of the precise nature of the statistics of 
the velocity field, the drift and diffusivity parameters are renormalized in the same 
way if we start from a situation where the corresponding microscopic quantities are 
proportional. This result confirms the theoretical prediction previously obtained in 
|I0| using a theoretical approach developed for the continuum from a method due 
to Derrida |], 111 O]. 



4 Conclusions 

In this paper we have studied the motion of a neutral molecule in a random gradient 
flow. The suggested physical mechanism giving rise to the local drift involves the 
interaction of the local electric field (presumed to have Gaussian statistics) with the 
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field-induced dipole moment of tlie molecule. The resulting drift field therefore does 
not have Gaussian statistics. That such a physically plausible model can give rise, 
in a natural way, to a drift field with non-Gaussian statistics, is itself interesting. 

The increased complexity of the model is no barrier to the formulation of a 
perturbative calculation of its Green's functions and effective parameters. The non- 
Gaussian character of the drift field statistics shows up first at three loop order. It 
results in a contribution to the effective diffusivity that is not symmetric under a 
change of sign of the coupling. This asymmetry is very clear in both the perturbation 
theory and numerical results exhibited in figure ^ . The same figure also exhibits 
the results of the numerical simulation. It shows that the inclusion of the four 
loop terms allows the perturbation series to give a reasonably good account of the 
effective diffusivity for the coupling parameter in the range —0.7 < Aq < 0.7 . 

The complexity of the non-Gaussian model, however, has so far prevented, a 
satisfactory formulation of the kinds of self-consistent perturbative calculation or 
renormalization group calculation that were rather successful for the Gaussian model 
P, 1^, 1^. In both cases the problem centres round the new types of vertices, 
not present in the original perturbation theory scheme, that are induced by the 
loop contributions of perturbation theory. The formulation of effective calculational 
schemes of this kind remains a goal of great interest since, as is clear from figure ^, 
low order perturbation theory is inadequate for situations of large disorder. Further 
investigations are in progress. 

In addition to the effective diffusivity the effective drift parameter, which con- 
trols the response of the molecules to an externally applied constant gradient field, 
is also a significant physical quantity. The success of the renormalization group 
calculation in the Gaussian model was in part due to the fact that it respected the 
Einstein relation, namely that the effective diffusivity and drift parameter are renor- 
malized by the same factor from their molecular values. The Einstein relation was 
also demonstrated in low order perturbation theory. In fact the Einstein relation 
was later shown to hold quite generally, independently of the statistical properties of 
the random medium, provided the molecular drift and diffusivity tensors were pro- 



portional to one another [10|: the Gaussian character of the model is not relelvant. 
It is therefore important that we have been able to give numerical confirmation (see 
table 1) of the validity of the Einstein relation between effective drift and diffusivity 
in the non-Gaussian model investigated in this paper. 

The introduction of directional effects into the drift and diffusivity as well as the 
statistics of the drift field, are further problems of great interest but considerably 
increased complexity. 
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Figure 1: Vertex diagram. 
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k-q-p 



Figure 2: Two-loop contribution to S(k). 
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Figure 3: Three-loop contribution to S(k). 
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k k-p-p' k-p-p'-q-q' k-p-q k 

(e) 

Figure 4: Four-loop contributions to S(k). 
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-1.5 -1 -0.5 0.5 1 1.5 

disorder strength 

Figure 5: Kg versus Aq assuming kq = 1. The simulation data arc siiown (o) to 
be compared with the prediction of two- loop (dot-dashed) , three- loop (dashed) and 
four-loop (solid) perturbation theory 
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